%% Michael D. Konig, Andrei Levchenko, Tim Rogers and Fabrizio Zilibotti
%% Aggregate fluctuations in adaptive production networks
%% Forthcoming in Proceedings of the National Academy of Sciences, 2022

%% Replication file (Matlab) for Figure 6, main text.

%% Requirements:

% "firm_factset_ids_sic_nation.csv" -> CSV file with the following colums: id; firm_factset_id; firm_sic; firm_nation
% "2015_adjacency_matrix.mat" -> Adjacency matrix (array) stored in Matlab data file (for the year 2015).
% "mle_eta.mat" -> Matlab file with estimated firm fixed effects. Outpt generated from "mle_attachment.m".
% "in_deg_distribution.mat" -> Matlab file with the empirical in-degree (number of suppliers) distribtion. 
% "European_countries.csv" -> CSV file with a list of European countriesand the following colums: Name,2 Letter Country Code,3 Letter Country Code,EU
% "graphconncomp.m" -> Find strongly or weakly connected components in graph. See: https://www.mathworks.com/help/bioinfo/ref/graphconncomp.html.
% "scomponents.m" -> Compute the strongly connected components of a graph. David F. Gleich, Copyright, Stanford University, 2008-2009.
% "charpath.m" -> Characteristic path length, global efficiency and related statistics.  Olaf Sporns, Indiana University, 2002/2007/2008 Modification, 2010 (Mika Rubinov): incorporation of global efficiency
% "modularity_dir.m" -> The optimal community structure is a subdivision of the network into non-overlapping groups of nodes in a way that maximizes the number of within-group edges, and minimizes the number of between-group edges. The modularity is a statistic that quantifies the degree to which the network may be subdivided into such clearly delineated groups. See: https://github.com/jblocher/matlab-network-utilities and http://strategic.mit.edu/downloads.php?page=matlab_networks.

clear all;
close all;
 
%% Parameters.
N = 300; % Input-Output characteristics / Binomial fitness distribution: number of trials.
emp_gamma = 0.1937; % Rewiring probability.
gamma = [[0 0.4 0.6 0.8 1.0 1.2 1.4] emp_gamma]; 
q = 0.375; % Output characteristics / Binomial fitness distribution: success probability.
p = 0.075; % Input characteristics success probability.
 
gamma_plot_min = min(gamma); 
gamma_plot_max = max(gamma);
 
phi = 0.6322; % Bonacich centrality parameter.
 
load_eta = true; % If true, load firm fixed effects (faitness values).
 
aggr_sic_level = 2;
drop_sic_missing = true;
 
dt = 0.0018; % Time increment.
T = 250; % Maximum number of iterations.
S = 500; % Number of realizations.
 
selected_year = 2015;
load_autarky = false;
 
multi_plot = true;
 
compute_connectivity_stats = true;
compute_strong_components = true; % Compute the strongly connected components. 
compute_weak_components = true; % Compute the weakly connected components.
compute_char_path_length = false; % Characteristic path length.
compute_modularity = false; % Modularity a la Newman and Leicht (The community structure is a subdivision of the network into nonoverlapping groups of nodes in a way that maximizes the number of within-group edges, and minimizes the number of between-group edges).
 
nearest_neighbors = true; % Average over the nearest neighbors when plotting the simulated sales coefficient of variation.
 
tic
 
rng(1,'twister')
 
%% Load firm id mappings data.
fid = fopen(['Data/firm_factset_ids_sic_nation.csv']);
dat = textscan(fid,repmat('%s ',1,4), 'delimiter', ';', 'headerlines', 1); % id; firm_factset_id; firm_sic; firm_nation
fclose(fid);
factset_firm_id = str2double(strtrim(dat{2}));
factset_firm_sic_str = strtrim(dat{3});
factset_firm_nation = strtrim(dat{4});
clear('dat')
factset_firm_nation_init = factset_firm_nation;
 
%% Aggregate SIC codes.
factset_firm_sic = nan(length(factset_firm_sic_str),1);
for i=1:length(factset_firm_sic_str)
    if(~strcmp(factset_firm_sic_str{i},'NaN'))
        dummy = factset_firm_sic_str{i};
        dummy = dummy(1:aggr_sic_level);
        if(~isempty(str2num(dummy)))
            factset_firm_sic(i) = str2num(dummy);
        end
    end
end
 
%% Load adjacency matrix.
disp(['Load adjacency matrix from file...'])
load (['./adjacency_matrix/' num2str(selected_year) '_adjacency_matrix.mat'])
 
%% Load firm fixed effects (fitness values).
if(load_eta)
    disp(['Load firm fixed effects (fitness values) from file...'])
    load (['./Data/mle_eta.mat']);
    eta = zeros(length(factset_firm_id),1);
    mle_eta_id = mle_eta(:,1);
    mle_eta_fit = mle_eta(:,2);
    for i=1:length(factset_firm_id)
        ind = find(factset_firm_id(i)==mle_eta_id);
        if(~isempty(ind))
            eta(i) = mle_eta_fit(ind(1));
        end
    end
    clear('mle_eta')
    mle_eta = [factset_firm_id eta];
else
    eta = zeros(length(factset_firm_id),1);
    for i=1:length(factset_firm_id)
        k = binornd(N,q);
        a = (1-p)^N;
        b = 1 + p/(q*(1-p));
        eta(i) = a*(b^k-1);
    end
    mle_eta = [factset_firm_id eta];
end
clear('eta')
 
%% Drop firms with missing SIC codes.
if(drop_sic_missing)
    disp(['Drop firms with missing SIC codes...'])
    ind = find(~isnan(factset_firm_sic));
    factset_firm_id = factset_firm_id(ind);
    factset_firm_sic = factset_firm_sic(ind);
    factset_firm_nation = factset_firm_nation(ind);
    A = A(ind,ind);
    mle_eta = mle_eta(ind,:);
end
 
A_benchmark = A;
 
%% Estimated rewiring rates.
rewiring_gamma = 0.1937*ones(length(rewiring_sic),1);
 
%% Estimated exit rates.
lifetime_delta = 0.097479*ones(length(lifetime_sic),1); 
lifetime_zeta = 0.0060066*ones(length(lifetime_sic),1); 
lifetime_rho = 0.0082788*ones(length(lifetime_sic),1); 
 
%% Load empirical in-degree distribution.
disp(['Load empirical in-degree distribution from file...'])
load('./Data/in_deg_distribution.mat');
ind = find(in_deg_distribution(:,1)==selected_year);
emp_indeg_distribution = in_deg_distribution(end,2:end);
 
%% Load list of European countries.
disp('Load firms locations..')
fid = fopen(['./Data/European_countries.csv']); %
dat = textscan(fid,repmat('%q ',1,4), 'delimiter', ',', 'headerlines', 1); %Name,2 Letter Country Code,3 Letter Country Code,EU
fclose(fid);
country_iso_2 = strtrim(dat{2});
country_eu_member = str2double(strtrim(dat{4}));
clear('dat')
ind = find(country_eu_member==0);
country_iso_2(ind) = []; % Keep only EU member states.
eu_iso_2 = country_iso_2;
 
%% Removal all external links in the US, Europe and Eastern Asia.
if(~load_autarky)
    
    A_autarky = A;
    ind_us = find(strcmp(factset_firm_nation,'US'));
    ind_not_us = setdiff([1:1:length(factset_firm_nation)],ind_us)'; % C = setdiff(A,B) returns the data in A that is not in B, with no repetitions. C is in sorted order.
    ind_eu = [];
    for j=1:length(eu_iso_2)
        ind_eu = [ind_eu; find(strcmp(factset_firm_nation,eu_iso_2{j}))];
    end
    ind_not_eu = setdiff([1:1:length(factset_firm_nation)],ind_eu)';
    ind_china = find(strcmp(factset_firm_nation,'CN'));
    ind_taiwan = find(strcmp(factset_firm_nation,'TW'));
    ind_japan = find(strcmp(factset_firm_nation,'JP'));
    ind_korea = find(strcmp(factset_firm_nation,'KR'));
    ind_singapore = find(strcmp(factset_firm_nation,'SG'));
    ind_hong_kong = find(strcmp(factset_firm_nation,'HK'));
    ind_asia = [ind_china; ind_taiwan; ind_japan; ind_korea; ind_singapore; ind_hong_kong];
    ind_not_asia = setdiff([1:1:length(factset_firm_nation)],ind_asia)';
    
    disp(['Block non-US links...'])
    num_iterations = length(ind_us);
    increment = floor(num_iterations/50);
    for i=1:length(ind_us)
        A_autarky(ind_us(i),ind_not_us) = 0;
        A_autarky(ind_not_us,ind_us(i)) = 0;
        %% Display progress of the for loop.
        if(mod(i,increment)==0)
            fprintf('Progress...%.0f%%\n',round(100*i/num_iterations));
        end
    end
    
    disp(['Block non-EU links...'])
    num_iterations = length(ind_eu);
    increment = floor(num_iterations/50);
    for i=1:length(ind_eu)
        A_autarky(ind_eu(i),ind_not_eu) = 0;
        A_autarky(ind_not_eu,ind_eu(i)) = 0;
        %% Display progress of the for loop.
        if(mod(i,increment)==0)
            fprintf('Progress...%.0f%%\n',round(100*i/num_iterations));
        end
    end
    
    disp(['Block non-Asian links...'])
    num_iterations = length(ind_asia);
    increment = floor(num_iterations/50);
    for i=1:length(ind_asia)
        A_autarky(ind_asia(i),ind_not_asia) = 0;
        A_autarky(ind_not_asia,ind_asia(i)) = 0;
        %% Display progress of the for loop.
        if(mod(i,increment)==0)
            fprintf('Progress...%.0f%%\n',round(100*i/num_iterations));
        end
    end
    
    save('./Data/A_autarky.mat','A_autarky')
    
else
    
    load('./Data/A_autarky.mat')
    
end
 
%% Iterate over different scenarios.
disp(['Iterate over different scenarios.'])
scenario = {'benchmark'; 'autarky'};
n_scenario = [];
n_scenario_us = [];
n_scenario_asia = [];
n_scenario_eu = [];
bonacich_mean_scenario = [];
bonacich_log_sum_std_scenario = [];
bonacich_coeff_var_scenario = [];
bonacich_mean_scenario_us = [];
bonacich_log_sum_std_scenario_us = [];
bonacich_coeff_var_scenario_us = [];
bonacich_mean_scenario_asia = [];
bonacich_log_sum_std_scenario_asia = [];
bonacich_coeff_var_scenario_asia = [];
bonacich_mean_scenario_eu = [];
bonacich_log_sum_std_scenario_eu = [];
bonacich_coeff_var_scenario_eu = [];
out_d_std_scenario = [];
deg_scenario = [];
out_d_std_scenario_us = [];
deg_scenario_us = [];
out_d_std_scenario_asia = [];
deg_scenario_asia = [];
out_d_std_scenario_eu = [];
deg_scenario_eu = [];
if(compute_connectivity_stats)
    if(compute_weak_components) num_weak_comp_scenario = []; end
    if(compute_weak_components) size_largest_weak_comp_scenario = []; end
    if(compute_strong_components) num_strong_comp_scenario = []; end
    if(compute_strong_components) size_largest_strong_comp_scenario = []; end
    if(compute_char_path_length) char_path_length_scenario = []; end
    if(compute_modularity) modularity_scenario = []; end
    
    if(compute_weak_components) num_weak_comp_scenario_us = []; end
    if(compute_weak_components) size_largest_weak_comp_scenario_us = []; end
    if(compute_strong_components) num_strong_comp_scenario_us = []; end
    if(compute_strong_components) size_largest_strong_comp_scenario_us = []; end
    if(compute_char_path_length) char_path_length_scenario_us = []; end
    if(compute_modularity) modularity_scenario_us = []; end
    
    if(compute_weak_components) num_weak_comp_scenario_asia = []; end
    if(compute_weak_components) size_largest_weak_comp_scenario_asia = []; end
    if(compute_strong_components) num_strong_comp_scenario_asia = []; end
    if(compute_strong_components) size_largest_strong_comp_scenario_asia = []; end
    if(compute_char_path_length) char_path_length_scenario_asia = []; end
    if(compute_modularity) modularity_scenario_asia = []; end
    
    if(compute_weak_components) num_weak_comp_scenario_eu = []; end
    if(compute_weak_components) size_largest_weak_comp_scenario_eu = []; end
    if(compute_strong_components) num_strong_comp_scenario_eu = []; end
    if(compute_strong_components) size_largest_strong_comp_scenario_eu = []; end
    if(compute_char_path_length) char_path_length_scenario_eu = []; end
    if(compute_modularity) modularity_scenario_eu = []; end
end
clear('A')
for c=1:length(scenario)
    
    if(strcmp(scenario(c),'benchmark'))
        
        fprintf('Compute benchmark... \n');
        A_0 = A_benchmark;
        
    else
        
        fprintf('Compute counterfactual... \n');
        A_0 = A_autarky;
        
    end
    n_0 = length(A_0);
    disp(['Number of firms ' num2str(n_0) '...']);
    
    %% Iterate over different values for the rewiring rate gamma.
    disp(['Iterate over different values for the rewiring rate gamma...'])
    n_save = [];
    n_save_us = [];
    n_save_asia = [];
    n_save_eu = [];
    n_std_save = [];
    n_coeff_var_save = [];
    deg_save = [];
    in_d_std_save = [];
    out_d_std_save = [];
    deg_save_us = [];
    in_d_std_save_us = [];
    out_d_std_save_us = [];
    deg_save_asia = [];
    in_d_std_save_asia = [];
    out_d_std_save_asia = [];
    deg_save_eu = [];
    in_d_std_save_eu = [];
    out_d_std_save_eu = [];
    if(compute_connectivity_stats)
        if(compute_weak_components)
            num_weak_comp_save = [];
            size_largest_weak_comp_save = [];
            num_weak_comp_save_us = [];
            size_largest_weak_comp_save_us = [];
            num_weak_comp_save_asia = [];
            size_largest_weak_comp_save_asia = [];
            num_weak_comp_save_eu = [];
            size_largest_weak_comp_save_eu = [];
        end
        if(compute_strong_components)
            num_strong_comp_save = [];
            size_largest_strong_comp_save = [];
            num_strong_comp_save_us = [];
            size_largest_strong_comp_save_us = [];
            num_strong_comp_save_asia = [];
            size_largest_strong_comp_save_asia = [];
            num_strong_comp_save_eu = [];
            size_largest_strong_comp_save_eu = [];
        end
        if(compute_char_path_length) char_path_length_save = []; end
        if(compute_modularity) modularity_save = []; end
        if(compute_char_path_length) char_path_length_save_us = []; end
        if(compute_modularity) modularity_save_us = []; end
        if(compute_char_path_length) char_path_length_save_asia = []; end
        if(compute_modularity) modularity_save_asia = []; end
        if(compute_char_path_length) char_path_length_save_eu = []; end
        if(compute_modularity) modularity_save_eu = []; end
    end
    bonacich_log_sum_save = [];
    bonacich_log_sum_std_save = [];
    bonacich_log_sum_skewness_save = [];
    bonacich_coeff_var_save = [];
    bonacich_mean_save = [];
    bonacich_std_save = [];
    bonacich_skewness_save = [];
    bonacich_std_growth_save = [];
    bonacich_log_sum_save_us = [];
    bonacich_log_sum_std_save_us = [];
    bonacich_log_sum_skewness_save_us = [];
    bonacich_coeff_var_save_us = [];
    bonacich_mean_save_us = [];
    bonacich_std_save_us = [];
    bonacich_skewness_save_us = [];
    bonacich_std_growth_save_us = [];
    bonacich_log_sum_save_asia = [];
    bonacich_log_sum_std_save_asia = [];
    bonacich_log_sum_skewness_save_asia = [];
    bonacich_coeff_var_save_asia = [];
    bonacich_mean_save_asia = [];
    bonacich_std_save_asia = [];
    bonacich_skewness_save_asia = [];
    bonacich_std_growth_save_asia = [];
    bonacich_log_sum_save_eu = [];
    bonacich_log_sum_std_save_eu = [];
    bonacich_log_sum_skewness_save_eu = [];
    bonacich_coeff_var_save_eu = [];
    bonacich_mean_save_eu = [];
    bonacich_std_save_eu = [];
    bonacich_skewness_save_eu = [];
    bonacich_std_growth_save_eu = [];
    
    for g=1:length(gamma)
        
        disp(['gamma = ' num2str(gamma(g))])
        
        %% Initialize vectors.
        n_time = zeros(T,S);
        n_time_us = zeros(T,S);
        n_time_asia = zeros(T,S);
        n_time_eu = zeros(T,S);
        deg_av_time = zeros(T,S);
        in_d_std_time = zeros(T,S);
        out_d_std_time = zeros(T,S);
        deg_av_time_us = zeros(T,S);
        in_d_std_time_us = zeros(T,S);
        out_d_std_time_us = zeros(T,S);
        deg_av_time_asia = zeros(T,S);
        in_d_std_time_asia = zeros(T,S);
        out_d_std_time_asia = zeros(T,S);
        deg_av_time_eu = zeros(T,S);
        in_d_std_time_eu = zeros(T,S);
        out_d_std_time_eu = zeros(T,S);
        bonacich_sum_time = zeros(T,S);
        bonacich_log_sum_time = zeros(T,S);
        bonacich_std_time = zeros(T,S);
        bonacich_skewness_time = zeros(T,S);
        bonacich_growth = zeros(T,S);
        bonacich_sum_time_us = zeros(T,S);
        bonacich_log_sum_time_us = zeros(T,S);
        bonacich_std_time_us = zeros(T,S);
        bonacich_skewness_time_us = zeros(T,S);
        bonacich_growth_us = zeros(T,S);
        bonacich_sum_time_asia = zeros(T,S);
        bonacich_log_sum_time_asia = zeros(T,S);
        bonacich_std_time_asia = zeros(T,S);
        bonacich_skewness_time_asia = zeros(T,S);
        bonacich_growth_asia = zeros(T,S);
        bonacich_sum_time_eu = zeros(T,S);
        bonacich_log_sum_time_eu = zeros(T,S);
        bonacich_std_time_eu = zeros(T,S);
        bonacich_skewness_time_eu = zeros(T,S);
        bonacich_growth_eu = zeros(T,S);
        if(compute_connectivity_stats)
            if(compute_weak_components) num_weak_comp_time = zeros(T,S); end
            if(compute_weak_components) size_largest_weak_comp_time = zeros(T,S); end
            if(compute_strong_components) num_strong_comp_time = zeros(T,S); end
            if(compute_strong_components) size_largest_strong_comp_time = zeros(T,S); end
            if(compute_char_path_length) char_path_length_time = zeros(T,S); end
            if(compute_modularity) modularity_time = zeros(T,S); end
            
            if(compute_weak_components) num_weak_comp_time_us = zeros(T,S); end
            if(compute_weak_components) size_largest_weak_comp_time_us = zeros(T,S); end
            if(compute_strong_components) num_strong_comp_time_us = zeros(T,S); end
            if(compute_strong_components) size_largest_strong_comp_time_us = zeros(T,S); end
            if(compute_char_path_length) char_path_length_time_us = zeros(T,S); end
            if(compute_modularity) modularity_time_us = zeros(T,S); end
            
            if(compute_weak_components) num_weak_comp_time_asia = zeros(T,S); end
            if(compute_weak_components) size_largest_weak_comp_time_asia = zeros(T,S); end
            if(compute_strong_components) num_strong_comp_time_asia = zeros(T,S); end
            if(compute_strong_components) size_largest_strong_comp_time_asia = zeros(T,S); end
            if(compute_char_path_length) char_path_length_time_asia = zeros(T,S); end
            if(compute_modularity) modularity_time_asia = zeros(T,S); end
            
            if(compute_weak_components) num_weak_comp_time_eu = zeros(T,S); end
            if(compute_weak_components) size_largest_weak_comp_time_eu = zeros(T,S); end
            if(compute_strong_components) num_strong_comp_time_eu = zeros(T,S); end
            if(compute_strong_components) size_largest_strong_comp_time_eu = zeros(T,S); end
            if(compute_char_path_length) char_path_length_time_eu = zeros(T,S); end
            if(compute_modularity) modularity_time_eu = zeros(T,S); end
        end
        
        %% Iterating over different realizations.
        for s = 1:S
            
            fprintf('Realization %3.0f.\n', s);
            
            %% Setting the seed for the random number generator.
            rand('seed',s);
            
            %% Generate initial network and assign initial fitness values.
            A = A_0;
            eta = mle_eta(:,2); % Assign initial fitness values.
            
            %% Assign delta, zeta and rho to incumbent firms.
            delta = 0.097479*ones(length(A_0),1); 
            zeta = 0.0060066*ones(length(A_0),1); 
            rho = 0.0082788*ones(length(A_0),1); 
            for i=1:length(A_0)
                ind = find(lifetime_sic==factset_firm_sic(i));
                if(~isempty(ind))
                    delta(i) = lifetime_delta(ind);
                    zeta(i) = lifetime_zeta(ind);
                    rho(i) = lifetime_rho(ind);
                end
            end
            
            %% Assign the number of initial suppliers.
            num_init_suppl = full(sum(A,1)); % in-degree.
            
            %% Assign US indicator.
            is_us = zeros(length(A),1);
            ind_us = find(strcmp(factset_firm_nation,'US'));
            is_us(ind_us) = 1;
            
            %% Assign Asia indicator.
            is_china = zeros(length(A),1);
            ind_china = find(strcmp(factset_firm_nation,'CN'));
            is_china(ind_china) = 1;
            is_taiwan = zeros(length(A),1);
            ind_taiwan = find(strcmp(factset_firm_nation,'TW'));
            is_taiwan(ind_taiwan) = 1;
            is_japan = zeros(length(A),1);
            ind_japan = find(strcmp(factset_firm_nation,'JP'));
            is_japan(ind_japan) = 1;
            is_korea = zeros(length(A),1);
            ind_korea = find(strcmp(factset_firm_nation,'KR'));
            is_korea(ind_korea) = 1;
            is_singapore = zeros(length(A),1);
            ind_singapore = find(strcmp(factset_firm_nation,'SG'));
            is_singapore(ind_singapore) = 1;
            is_hong_kong = zeros(length(A),1);
            ind_hong_kong = find(strcmp(factset_firm_nation,'HK'));
            is_hong_kong(ind_hong_kong) = 1;
            is_asia = zeros(length(A),1);
            is_asia(ind_china) = 1;
            is_asia(ind_taiwan) = 1;
            is_asia(ind_japan) = 1;
            is_asia(ind_korea) = 1;
            is_asia(ind_singapore) = 1;
            is_asia(ind_hong_kong) = 1;
            
            %% Assign EU indicator.
            ind_eu = [];
            for j=1:length(eu_iso_2)
                ind_eu = [ind_eu; find(strcmp(factset_firm_nation,eu_iso_2{j}))];
            end
            is_eu = zeros(length(A),1);
            is_eu(ind_eu) = 1;
            
            %% Network formation.
            for t=1:T
                
                %% Compute the time varying statistics.
                n = length(A);
                u = ones(n,1);
                out_d = full(sum(A,2)); % out-degree.
                in_d = full(sum(A,1))'; % in-degree.
                deg_av_time(t,s) = mean(in_d);
                in_d_std_time(t,s) = std(in_d);
                out_d_std_time(t,s) = std(out_d);
                n_time(t,s) = length(A);
                d_inv = 1./in_d;
                d_inv(isinf(d_inv)) = 1;
                W = A.*d_inv;
                bonacich = u;
                bonacich = (eye(n) - phi*W)\u; % Replace inv(A)*b with A\b.
                bonacich_sum_time(t,s) = nansum(bonacich);
                bonacich_log_sum_time(t,s) = log(sum(bonacich)); % log(sum(bonacich)) sum(log(bonacich))
                bonacich_std_time(t,s) = std(bonacich);
                bonacich_skewness_time(t,s) = skewness(bonacich);
                if(t>1 && bonacich_log_sum_time(t,s)>0)
                    bonacich_growth(t-1,s) = (bonacich_log_sum_time(t,s)-bonacich_log_sum_time(t-1,s))/bonacich_log_sum_time(t-1,s);
                end
                if(compute_connectivity_stats)
                    if(compute_weak_components) % Compute the weakly connected components
                        [dummy1, dummy2] = graphconncomp(sparse(A),'Directed',1,'Weak',1);
                        num_weak_comp_time(t,s) = dummy1;
                        ind_unique = unique(dummy2);
                        count = countmember(ind_unique,dummy2);
                        [max_size max_ind] = max(count);
                        size_largest_weak_comp_time(t,s) = max_size/length(A);
                    end
                    if(compute_strong_components) % Compute the strongly connected components
                        [ci sizes] = scomponents(A);
                        [csize cind] = max(sizes);
                        dummy = ci==cind;
                        Ascc = A(dummy,dummy);
                        size_largest_strong_comp_time(t,s) = length(Ascc)/length(A);
                        %[ci sizes] = scomponents(A);
                        if(length(ci)>0)
                            num_strong_comp_time(t,s) = max(ci);
                        else
                            num_strong_comp_time(t,s) = 1;
                        end
                    end
                    if(compute_char_path_length) % Characteristic path length.
                        [char_path_length_time(t,s),efficiency,ecc,radius,diameter] = charpath(distance_bin(A));
                    end
                    if(compute_modularity) % Modularity.
                         [Ci Q] = modularity_dir(A);
                         modularity_time(t,s) = Q;
                    end
                end
                ind = find(is_us);
                n_time_us(t,s) = length(ind);
                bonacich_sum_time_us(t,s) = nansum(bonacich(ind));
                bonacich_log_sum_time_us(t,s) = log(sum(bonacich(ind))); % log(sum(bonacich)) sum(log(bonacich))
                bonacich_std_time_us(t,s) = std(bonacich(ind));
                bonacich_skewness_time_us(t,s) = skewness(bonacich(ind));
                if(t>1 && bonacich_log_sum_time_us(t,s)>0)
                    bonacich_growth_us(t-1,s) = (bonacich_log_sum_time_us(t,s)-bonacich_log_sum_time_us(t-1,s))/bonacich_log_sum_time_us(t-1,s);
                end
                deg_av_time_us(t,s) = mean(in_d(ind));
                in_d_std_time_us(t,s) = std(in_d(ind));
                out_d_std_time_us(t,s) = std(out_d(ind));
                if(compute_connectivity_stats)
                    A_us = A(ind,ind);
                    if(compute_weak_components) % Compute the weakly connected components
                        [dummy1, dummy2] = graphconncomp(sparse(A_us),'Directed',1,'Weak',1);
                        num_weak_comp_time_us(t,s) = dummy1;
                        ind_unique = unique(dummy2);
                        count = countmember(ind_unique,dummy2);
                        [max_size max_ind] = max(count);
                        size_largest_weak_comp_time_us(t,s) = max_size/length(A_us);
                    end
                    if(compute_strong_components) % Compute the strongly connected components
                        [ci sizes] = scomponents(A_us);
                        [csize cind] = max(sizes);
                        dummy = ci==cind;
                        Ascc_us = A_us(dummy,dummy);
                        size_largest_strong_comp_time_us(t,s) = length(Ascc_us)/length(A_us);
                        %[ci sizes] = scomponents(A_us);
                        if(length(ci)>0)
                            num_strong_comp_time_us(t,s) = max(ci);
                        else
                            num_strong_comp_time_us(t,s) = 1;
                        end
                    end
                    if(compute_char_path_length) % Characteristic path length.
                        [char_path_length_time_us(t,s),efficiency,ecc,radius,diameter] = charpath(distance_bin(A_us));
                    end
                    if(compute_modularity) % Modularity.
                         [Ci Q] = modularity_dir(A_us);
                         modularity_time_us(t,s) = Q;
                    end
                end
                ind = find(is_asia);
                n_time_asia(t,s) = length(ind);
                bonacich_sum_time_asia(t,s) = nansum(bonacich(ind));
                bonacich_log_sum_time_asia(t,s) = log(sum(bonacich(ind))); 
                bonacich_std_time_asia(t,s) = std(bonacich(ind));
                bonacich_skewness_time_asia(t,s) = skewness(bonacich(ind));
                if(t>1 && bonacich_log_sum_time_asia(t,s)>0)
                    bonacich_growth_asia(t-1,s) = (bonacich_log_sum_time_asia(t,s)-bonacich_log_sum_time_asia(t-1,s))/bonacich_log_sum_time_asia(t-1,s);
                end
                deg_av_time_asia(t,s) = mean(in_d(ind));
                in_d_std_time_asia(t,s) = std(in_d(ind));
                out_d_std_time_asia(t,s) = std(out_d(ind));
                if(compute_connectivity_stats)
                    A_asia = A(ind,ind);
                    if(compute_weak_components) % Compute the weakly connected components
                        [dummy1, dummy2] = graphconncomp(sparse(A_asia),'Directed',1,'Weak',1);
                        num_weak_comp_time_asia(t,s) = dummy1;
                        ind_unique = unique(dummy2);
                        count = countmember(ind_unique,dummy2);
                        [max_size max_ind] = max(count);
                        size_largest_weak_comp_time_asia(t,s) = max_size/length(A_asia);
                    end
                    if(compute_strong_components) % Compute the strongly connected components
                        [ci sizes] = scomponents(A_asia);
                        [csize cind] = max(sizes);
                        dummy = ci==cind;
                        Ascc_asia = A_asia(dummy,dummy);
                        size_largest_strong_comp_time_asia(t,s) = length(Ascc_asia)/length(A_asia);
                        %[ci sizes] = scomponents(A_asia);
                        if(length(ci)>0)
                            num_strong_comp_time_asia(t,s) = max(ci);
                        else
                            num_strong_comp_time_asia(t,s) = 1;
                        end
                    end
                    if(compute_char_path_length) % Characteristic path length.
                        [char_path_length_time_asia(t,s),efficiency,ecc,radius,diameter] = charpath(distance_bin(A_asia));
                    end
                    if(compute_modularity) % Modularity.
                         [Ci Q] = modularity_dir(A_asia);
                         modularity_time_asia(t,s) = Q;
                    end
                end
                ind = find(is_eu);
                n_time_eu(t,s) = length(ind);
                bonacich_sum_time_eu(t,s) = nansum(bonacich(ind));
                bonacich_log_sum_time_eu(t,s) = log(sum(bonacich(ind))); 
                bonacich_std_time_eu(t,s) = std(bonacich(ind));
                bonacich_skewness_time_eu(t,s) = skewness(bonacich(ind));
                if(t>1 && bonacich_log_sum_time_eu(t,s)>0)
                    bonacich_growth_eu(t-1,s) = (bonacich_log_sum_time_eu(t,s)-bonacich_log_sum_time_eu(t-1,s))/bonacich_log_sum_time_eu(t-1,s);
                end
                deg_av_time_eu(t,s) = mean(in_d(ind));
                in_d_std_time_eu(t,s) = std(in_d(ind));
                out_d_std_time_eu(t,s) = std(out_d(ind));
                if(compute_connectivity_stats)
                    A_eu = A(ind,ind);
                    if(compute_weak_components) % Compute the weakly connected components
                        [dummy1, dummy2] = graphconncomp(sparse(A_eu),'Directed',1,'Weak',1);
                        num_weak_comp_time_eu(t,s) = dummy1;
                        ind_unique = unique(dummy2);
                        count = countmember(ind_unique,dummy2);
                        [max_size max_ind] = max(count);
                        size_largest_weak_comp_time_eu(t,s) = max_size/length(A_eu);
                    end
                    if(compute_strong_components) % Compute the strongly connected components
                        [ci sizes] = scomponents(A_eu);
                        [csize cind] = max(sizes);
                        dummy = ci==cind;
                        Ascc_eu = A_eu(dummy,dummy);
                        size_largest_strong_comp_time_eu(t,s) = length(Ascc_eu)/length(A_eu);
                        %[ci sizes] = scomponents(A_eu);
                        if(length(ci)>0)
                            num_strong_comp_time_eu(t,s) = max(ci);
                        else
                            num_strong_comp_time_eu(t,s) = 1;
                        end
                    end
                    if(compute_char_path_length) % Characteristic path length.
                        [char_path_length_time_eu(t,s),efficiency,ecc,radius,diameter] = charpath(distance_bin(A_eu));
                    end
                    if(compute_modularity) % Modularity.
                         [Ci Q] = modularity_dir(A_eu);
                         modularity_time_eu(t,s) = Q;
                    end
                end
                
                %% With probability delta*dt, remove each node at random.
                if(length(A)>1)
                    ind_remove = [];
                    for k=1:length(A)
                        if(rand<delta(k)*dt)
                            ind_remove = [ind_remove k];
                        end
                    end
                    while(~isempty(ind_remove))
                        k = ind_remove(1);
                        ind_remove(1) = [];
                        ind = find(ind_remove>k);
                        for j=1:length(ind)
                            ind_remove(ind(j)) = ind_remove(ind(j))-1;
                        end
                        A(k,:) = [];
                        A(:,k) = [];
                        eta(k) = [];
                        delta(k) = [];
                        zeta(k) = [];
                        rho(k) = [];
                        num_init_suppl(k) = [];
                        is_us(k) = [];
                        is_asia(k) = [];
                        is_eu(k) = [];
                    end
                end
                
                %% With probability zeta*dt, remove each node at random if its after shock profits are negative (and profits are proportional to the out-degree).
                if(length(A)>1)
                    ind_remove = [];
                    for k=1:length(A)
                        %if(rand<zeta(k)*dt/(out_d(k)+1))
                        if(rand<zeta(k)*dt/bonacich(k))
                            ind_remove = [ind_remove k];
                        end
                    end
                    while(~isempty(ind_remove))
                        k = ind_remove(1);
                        ind_remove(1) = [];
                        ind = find(ind_remove>k);
                        for j=1:length(ind)
                            ind_remove(ind(j)) = ind_remove(ind(j))-1;
                        end
                        A(k,:) = [];
                        A(:,k) = [];
                        eta(k) = [];
                        delta(k) = [];
                        zeta(k) = [];
                        rho(k) = [];
                        num_init_suppl(k) = [];
                        is_us(k) = [];
                        is_asia(k) = [];
                        is_eu(k) = [];
                    end
                end
                
                %% With probability rho*dt select each firm which has an in-degree zero and remove it from the network.
                if(length(A)>1)
                    ind_remove = [];
                    for k=1:length(A)
                        if(in_d(k)==0 && rand<rho(k)*dt)
                            ind_remove = [ind_remove k];
                        end
                    end
                    while(~isempty(ind_remove))
                        k = ind_remove(1);
                        ind_remove(1) = [];
                        ind = find(ind_remove>k);
                        for j=1:length(ind)
                            ind_remove(ind(j)) = ind_remove(ind(j))-1;
                        end
                        A(k,:) = [];
                        A(:,k) = [];
                        eta(k) = [];
                        delta(k) = [];
                        zeta(k) = [];
                        rho(k) = [];
                        num_init_suppl(k) = [];
                        is_us(k) = [];
                        is_asia(k) = [];
                        is_eu(k) = [];
                        for l=1:length(ind_remove)
                            if(ind_remove(l)>k)
                                ind_remove(l) = ind_remove(l)-1;
                            end
                        end
                    end
                end
                
                %% With probability gamma*dt select each firm which has in-degree less than the inital in-degree and find a new supplier from the incumbent firms (rewiring).
                if(length(A)>1)
                    ind_rewire = [];
                    for k=1:length(A)
                        if(in_d(k)< num_init_suppl(k) && rand<gamma(g)*dt) %if(in_d(k)==0 && rand<gamma(k)*dt)
                            ind_rewire = [ind_rewire k];
                        end
                    end
                    n = length(A);
                    while(~isempty(ind_rewire))
                        
                        ind = ind_rewire(1);
                        ind_rewire(1) = [];
                        
                        %% Fitness (eta) based attachment.
                        if(strcmp(scenario(c),'autarky'))
                            if(is_us(ind))
                                dummy = find(is_us==1);
                            elseif(is_eu(ind))
                                dummy = find(is_eu==1);
                            elseif(is_asia(ind))
                                dummy = find(is_asia==1);
                            else
                                dummy = [find(is_us==1); find(is_eu==1); find(is_asia==1)];
                                dummy = setdiff([1:1:n],dummy)'; % C = setdiff(A,B) returns the data in A that is not in B, with no repetitions. C is in sorted order.
                            end
                            dummy(dummy==ind) = [];
                            j = 1;
                            f = eta(dummy)/length(dummy);
                            f = f/sum(f);
                            F = f(j);
                            u = rand;
                            while(F<u)
                                j = j+1;
                                F = F + f(j);
                            end
                            A(dummy(j),ind) = 1;
                        else
                            j = 1;
                            f = eta/n;
                            f = f/sum(f);
                            F = f(j);
                            u = rand;
                            while(F<u)
                                j = j+1;
                                F = F + f(j);
                            end
                            A(j,ind) = 1;
                        end
                        
                    end
                end
                
                %% Add a new node.
                n = length(A);
                A = [A; sparse(1,n)]; % Expanding the adjacency matrix by one column and row of zeros.
                A = [A, sparse(n+1,1)];
                n = length(A);
                
                %% Assign new fitness value.
                k = binornd(N,q);
                a = (1-p)^N;
                b = 1 + p/(q*(1-p));
                eta(n) = a*(b^k-1);
                
                %% Assign rates.
                delta(n) = lifetime_delta(randi(length(lifetime_delta)));
                zeta(n) = lifetime_zeta(randi(length(lifetime_zeta)));
                rho(n) = lifetime_rho(randi(length(lifetime_rho)));
                
                %% Assign country.
                %country = factset_firm_nation(randi(length(factset_firm_nation)));
                country = factset_firm_nation_init(randi(length(factset_firm_nation_init)));
                if(strcmp(country,'US'))
                    is_us = [is_us; 1];
                    is_asia = [is_asia; 0];
                    is_eu = [is_eu; 0];
                    ind_us = [ind_us; n];
                elseif(strcmp(country,'CN'))
                    is_us = [is_us; 0];
                    is_eu = [is_eu; 0];
                    is_asia = [is_asia; 1];
                    ind_asia = [ind_asia; n];
                elseif(strcmp(country,'TW'))
                    is_us = [is_us; 0];
                    is_eu = [is_eu; 0];
                    is_asia = [is_asia; 1];
                    ind_asia = [ind_asia; n];
                elseif(strcmp(country,'JP'))
                    is_us = [is_us; 0];
                    is_eu = [is_eu; 0];
                    is_asia = [is_asia; 1];
                    ind_asia = [ind_asia; n];                    
                elseif(strcmp(country,'KR'))
                    is_us = [is_us; 0];
                    is_eu = [is_eu; 0];
                    is_asia = [is_asia; 1];
                    ind_asia = [ind_asia; n];                     
                elseif(strcmp(country,'SG'))
                    is_us = [is_us; 0];
                    is_eu = [is_eu; 0];
                    is_asia = [is_asia; 1];
                    ind_asia = [ind_asia; n];                    
                elseif(strcmp(country,'HK'))
                    is_us = [is_us; 0];
                    is_eu = [is_eu; 0];
                    is_asia = [is_asia; 1];
                    ind_asia = [ind_asia; n]; 
                else
                    found_eu = 0;
                    for j=1:length(eu_iso_2)
                        if(strcmp(country,eu_iso_2{j}))
                            is_us = [is_us; 0];
                            is_asia = [is_asia; 0];
                            is_eu = [is_eu; 1];
                            ind_eu = [ind_eu; n];
                            found_eu = 1;
                        end
                    end
                    if(~found_eu)
                        is_us = [is_us; 0];
                        is_asia = [is_asia; 0];
                        is_eu = [is_eu; 0];
                    end
                end
                
                %% Draw the number of suppliers of the firm.
                l = 1;
                F = emp_indeg_distribution(l);
                u = rand;
                while(F<u)
                    l = l+1;
                    F = F + emp_indeg_distribution(l);
                end
                num_init_suppl(n) = l; % Record the number of initial suppliers.
                
                %% Draw l suppliers from the incumbent firms.
                for k=1:num_init_suppl(n)
                    
                    %% Fitness (eta) based attachment.
                    if(strcmp(scenario(c),'autarky'))
                        if(is_us(n))
                            dummy = find(is_us==1);
                        elseif(is_eu(n))
                            dummy = find(is_eu==1);
                        elseif(is_asia(n))
                            dummy = find(is_asia==1);
                         else
                            dummy = [find(is_us==1); find(is_eu==1); find(is_asia==1)];
                            dummy = setdiff([1:1:n],dummy)'; % C = setdiff(A,B) returns the data in A that is not in B, with no repetitions. C is in sorted order.
                        end
                        dummy(dummy==n) = [];
                        j = 1;
                        f = eta(dummy)/length(dummy);
                        F = f(j);
                        u = rand;
                        while(F<u)
                            j = j+1;
                            F = F + f(j);
                        end
                        A(dummy(j),n) = 1;
                    else
                        j = 1;
                        f = eta/n;
                        f = f/sum(f);
                        F = f(j);
                        u = rand;
                        while(F<u)
                            j = j+1;
                            F = F + f(j);
                        end
                        A(j,n) = 1;
                    end
                    
                end
                
            end
            
        end
        
        n_time_mean = nanmean(n_time(floor(T/2):end,:),2);
        n_save = [n_save, nanmean(n_time_mean)];
        n_std = std(n_time(floor(T/2):end,:),0,1);
        n_std_save = [n_std_save, nanmean(n_std)];
        n_coeff_var_save = [n_coeff_var_save, nanmean(n_std)/mean(n_time_mean)];
        
        n_time_mean_us = nanmean(n_time_us(floor(T/2):end,:),2);
        n_save_us = [n_save_us, nanmean(n_time_mean_us)];
        
        n_time_mean_asia = nanmean(n_time_asia(floor(T/2):end,:),2);
        n_save_asia = [n_save_asia, nanmean(n_time_mean_asia)];
        
        n_time_mean_eu = nanmean(n_time_eu(floor(T/2):end,:),2);
        n_save_eu = [n_save_eu, nanmean(n_time_mean_eu)];
        
        deg_time_mean = mean(deg_av_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        deg_save = [deg_save, mean(deg_time_mean)]; % Average across time t=1,...,T.
        in_d_std_time_mean = mean(in_d_std_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        in_d_std_save = [in_d_std_save, mean(in_d_std_time_mean)]; % Average across time t=1,...,T.
        out_d_std_time_mean = mean(out_d_std_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        out_d_std_save = [out_d_std_save, nanmean(out_d_std_time_mean)]; % Average across time t=1,...,T.
        
        deg_time_mean_us = mean(deg_av_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        deg_save_us = [deg_save_us, mean(deg_time_mean_us)]; % Average across time t=1,...,T.
        in_d_std_time_mean_us = mean(in_d_std_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        in_d_std_save_us = [in_d_std_save_us, mean(in_d_std_time_mean_us)]; % Average across time t=1,...,T.
        out_d_std_time_mean_us = mean(out_d_std_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        out_d_std_save_us = [out_d_std_save_us, nanmean(out_d_std_time_mean_us)]; % Average across time t=1,...,T.
        
        if(compute_connectivity_stats)
            if(compute_weak_components)
                num_weak_comp_time_mean = mean(num_weak_comp_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_weak_comp_save = [num_weak_comp_save, mean(num_weak_comp_time_mean)]; % Average across time t=1,...,T.
                size_largest_weak_comp_time_mean = mean(size_largest_weak_comp_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_weak_comp_save = [size_largest_weak_comp_save, mean(size_largest_weak_comp_time_mean)]; % Average across time t=1,...,T.
                
                num_weak_comp_time_mean_us = mean(num_weak_comp_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_weak_comp_save_us = [num_weak_comp_save_us, mean(num_weak_comp_time_mean_us)]; % Average across time t=1,...,T.
                size_largest_weak_comp_time_mean_us = mean(size_largest_weak_comp_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_weak_comp_save_us = [size_largest_weak_comp_save_us, mean(size_largest_weak_comp_time_mean_us)]; % Average across time t=1,...,T.
                
                num_weak_comp_time_mean_asia = mean(num_weak_comp_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_weak_comp_save_asia = [num_weak_comp_save_asia, mean(num_weak_comp_time_mean_asia)]; % Average across time t=1,...,T.
                size_largest_weak_comp_time_mean_asia = mean(size_largest_weak_comp_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_weak_comp_save_asia = [size_largest_weak_comp_save_asia, mean(size_largest_weak_comp_time_mean_asia)]; % Average across time t=1,...,T.
                
                num_weak_comp_time_mean_eu = mean(num_weak_comp_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_weak_comp_save_eu = [num_weak_comp_save_eu, mean(num_weak_comp_time_mean_eu)]; % Average across time t=1,...,T.
                size_largest_weak_comp_time_mean_eu = mean(size_largest_weak_comp_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_weak_comp_save_eu = [size_largest_weak_comp_save_eu, mean(size_largest_weak_comp_time_mean_eu)]; % Average across time t=1,...,T.
            end
            if(compute_strong_components)
                num_strong_comp_time_mean = mean(num_strong_comp_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_strong_comp_save = [num_strong_comp_save, mean(num_strong_comp_time_mean)]; % Average across time t=1,...,T.
                size_largest_strong_comp_time_mean = mean(size_largest_strong_comp_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_strong_comp_save = [size_largest_strong_comp_save, mean(size_largest_strong_comp_time_mean)]; % Average across time t=1,...,T.
 
                num_strong_comp_time_mean_us = mean(num_strong_comp_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_strong_comp_save_us = [num_strong_comp_save_us, mean(num_strong_comp_time_mean_us)]; % Average across time t=1,...,T.
                size_largest_strong_comp_time_mean_us = mean(size_largest_strong_comp_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_strong_comp_save_us = [size_largest_strong_comp_save_us, mean(size_largest_strong_comp_time_mean_us)]; % Average across time t=1,...,T.
 
                num_strong_comp_time_mean_asia = mean(num_strong_comp_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_strong_comp_save_asia = [num_strong_comp_save_asia, mean(num_strong_comp_time_mean_asia)]; % Average across time t=1,...,T.
                size_largest_strong_comp_time_mean_asia = mean(size_largest_strong_comp_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_strong_comp_save_asia = [size_largest_strong_comp_save_asia, mean(size_largest_strong_comp_time_mean_asia)]; % Average across time t=1,...,T.
 
                num_strong_comp_time_mean_eu = mean(num_strong_comp_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                num_strong_comp_save_eu = [num_strong_comp_save_eu, mean(num_strong_comp_time_mean_eu)]; % Average across time t=1,...,T.
                size_largest_strong_comp_time_mean_eu = mean(size_largest_strong_comp_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                size_largest_strong_comp_save_eu = [size_largest_strong_comp_save_eu, mean(size_largest_strong_comp_time_mean_eu)]; % Average across time t=1,...,T.
                
            end
            if(compute_char_path_length)
                char_path_length_time_mean = mean(char_path_length_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                char_path_length_save = [char_path_length_save, mean(char_path_length_time_mean)]; % Average across time t=1,...,T.
                
                char_path_length_time_mean_us = mean(char_path_length_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                char_path_length_save_us = [char_path_length_save_us, mean(char_path_length_time_mean_us)]; % Average across time t=1,...,T.
                
                char_path_length_time_mean_asia = mean(char_path_length_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                char_path_length_save_asia = [char_path_length_save_asia, mean(char_path_length_time_mean_asia)]; % Average across time t=1,...,T.
                
                char_path_length_time_mean_eu = mean(char_path_length_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                char_path_length_save_eu = [char_path_length_save_eu, mean(char_path_length_time_mean_eu)]; % Average across time t=1,...,T.
            end
            if(compute_modularity)
                modularity_time_mean = mean(modularity_time(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                modularity_save = [modularity_save, mean(modularity_time_mean)]; % Average across time t=1,...,T.
                
                modularity_time_mean_us = mean(modularity_time_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                modularity_save_us = [modularity_save_us, mean(modularity_time_mean_us)]; % Average across time t=1,...,T.
                
                modularity_time_mean_asia = mean(modularity_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                modularity_save_asia = [modularity_save_asia, mean(modularity_time_mean_asia)]; % Average across time t=1,...,T.
                
                modularity_time_mean_eu = mean(modularity_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
                modularity_save_eu = [modularity_save_eu, mean(modularity_time_mean_eu)]; % Average across time t=1,...,T.
            end
        end
        
        deg_time_mean_asia = mean(deg_av_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        deg_save_asia = [deg_save_asia, mean(deg_time_mean_asia)]; % Average across time t=1,...,T.
        in_d_std_time_mean_asia = mean(in_d_std_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        in_d_std_save_asia = [in_d_std_save_asia, mean(in_d_std_time_mean_asia)]; % Average across time t=1,...,T.
        out_d_std_time_mean_asia = mean(out_d_std_time_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        out_d_std_save_asia = [out_d_std_save_asia, nanmean(out_d_std_time_mean_asia)]; % Average across time t=1,...,T.
 
        deg_time_mean_eu = mean(deg_av_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        deg_save_eu = [deg_save_eu, mean(deg_time_mean_eu)]; % Average across time t=1,...,T.
        in_d_std_time_mean_eu = mean(in_d_std_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        in_d_std_save_eu = [in_d_std_save_eu, mean(in_d_std_time_mean_eu)]; % Average across time t=1,...,T.
        out_d_std_time_mean_eu = mean(out_d_std_time_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        out_d_std_save_eu = [out_d_std_save_eu, nanmean(out_d_std_time_mean_eu)]; % Average across time t=1,...,T.
 
        bonacich_log_sum_time_mean = nanmean(bonacich_log_sum_time(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_log_sum_save = [bonacich_log_sum_save, nanmean(bonacich_log_sum_time_mean)]; % Average across simulation runs s=1,...,S.
        bonacich_log_sum_std = nanstd(bonacich_log_sum_time(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_log_sum_std_save = [bonacich_log_sum_std_save, nanmean(bonacich_log_sum_std)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_log_sum_skewness = skewness(bonacich_log_sum_time(floor(T/2):end,:),1,1); % Skewness across time t=1,...,T.
        bonacich_log_sum_skewness_save = [bonacich_log_sum_skewness_save, nanmean(bonacich_log_sum_skewness)]; % Average acrosssimulation runs s=1,...,S.
        
        bonacich_mean = nanmean(bonacich_sum_time(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_mean_save = [bonacich_mean_save, nanmean(bonacich_mean)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_std = nanstd(bonacich_sum_time(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_save = [bonacich_std_save, nanmean(bonacich_std)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_coeff_var_save = [bonacich_coeff_var_save, nanmean(bonacich_std./bonacich_mean)]; % Average acrosssimulation runs s=1,...,S.
        
        %bonacich_mean_growth = nanmean(bonacich_growth(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        bonacich_std_growth = std(bonacich_growth(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_growth_save = [bonacich_std_growth_save, nanmean(bonacich_std_growth)]; % Average across simulation runs s=1,...,S.
        
        bonacich_log_sum_time_mean_us = nanmean(bonacich_log_sum_time_us(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_log_sum_save_us = [bonacich_log_sum_save_us, nanmean(bonacich_log_sum_time_mean_us)]; % Average across simulation runs s=1,...,S.
        bonacich_log_sum_std_us = nanstd(bonacich_log_sum_time_us(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_log_sum_std_save_us = [bonacich_log_sum_std_save_us, nanmean(bonacich_log_sum_std_us)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_log_sum_skewness_us = skewness(bonacich_log_sum_time_us(floor(T/2):end,:),1,1); % Skewness across time t=1,...,T.
        bonacich_log_sum_skewness_save_us = [bonacich_log_sum_skewness_save_us, nanmean(bonacich_log_sum_skewness_us)]; % Average acrosssimulation runs s=1,...,S.
        
        bonacich_mean_us = nanmean(bonacich_sum_time_us(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_mean_save_us = [bonacich_mean_save_us, nanmean(bonacich_mean_us)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_std_us = nanstd(bonacich_sum_time_us(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_save_us = [bonacich_std_save_us, nanmean(bonacich_std_us)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_coeff_var_save_us = [bonacich_coeff_var_save_us, nanmean(bonacich_std_us./bonacich_mean_us)]; % Average acrosssimulation runs s=1,...,S.
        
        %bonacich_mean_growth_us = nanmean(bonacich_growth_us(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        bonacich_std_growth_us = std(bonacich_growth_us(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_growth_save_us = [bonacich_std_growth_save_us, nanmean(bonacich_std_growth_us)]; % Average across simulation runs s=1,...,S.
        
        bonacich_log_sum_time_mean_asia = nanmean(bonacich_log_sum_time_asia(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_log_sum_save_asia = [bonacich_log_sum_save_asia, nanmean(bonacich_log_sum_time_mean_asia)]; % Average across simulation runs s=1,...,S.
        bonacich_log_sum_std_asia = nanstd(bonacich_log_sum_time_asia(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_log_sum_std_save_asia = [bonacich_log_sum_std_save_asia, nanmean(bonacich_log_sum_std_asia)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_log_sum_skewness_asia = skewness(bonacich_log_sum_time_asia(floor(T/2):end,:),1,1); % Skewness across time t=1,...,T.
        bonacich_log_sum_skewness_save_asia = [bonacich_log_sum_skewness_save_asia, nanmean(bonacich_log_sum_skewness_asia)]; % Average acrosssimulation runs s=1,...,S.
        
        bonacich_mean_asia = nanmean(bonacich_sum_time_asia(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_mean_save_asia = [bonacich_mean_save_asia, nanmean(bonacich_mean_asia)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_std_asia = nanstd(bonacich_sum_time_asia(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_save_asia = [bonacich_std_save_asia, nanmean(bonacich_std_asia)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_coeff_var_save_asia = [bonacich_coeff_var_save_asia, nanmean(bonacich_std_asia./bonacich_mean_asia)]; % Average acrosssimulation runs s=1,...,S.
        
        %bonacich_mean_growth_asia = nanmean(bonacich_growth_asia(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        bonacich_std_growth_asia = std(bonacich_growth_asia(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_growth_save_asia = [bonacich_std_growth_save_asia, nanmean(bonacich_std_growth_asia)]; % Average across simulation runs s=1,...,S.
 
        bonacich_log_sum_time_mean_eu = nanmean(bonacich_log_sum_time_eu(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_log_sum_save_eu = [bonacich_log_sum_save_eu, nanmean(bonacich_log_sum_time_mean_eu)]; % Average across simulation runs s=1,...,S.
        bonacich_log_sum_std_eu = nanstd(bonacich_log_sum_time_eu(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_log_sum_std_save_eu = [bonacich_log_sum_std_save_eu, nanmean(bonacich_log_sum_std_eu)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_log_sum_skewness_eu = skewness(bonacich_log_sum_time_eu(floor(T/2):end,:),1,1); % Skewness across time t=1,...,T.
        bonacich_log_sum_skewness_save_eu = [bonacich_log_sum_skewness_save_eu, nanmean(bonacich_log_sum_skewness_eu)]; % Average acrosssimulation runs s=1,...,S.
        
        bonacich_mean_eu = nanmean(bonacich_sum_time_eu(floor(T/2):end,:),1); % Average across time t=1,...,T.
        bonacich_mean_save_eu = [bonacich_mean_save_eu, nanmean(bonacich_mean_eu)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_std_eu = nanstd(bonacich_sum_time_eu(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_save_eu = [bonacich_std_save_eu, nanmean(bonacich_std_eu)]; % Average acrosssimulation runs s=1,...,S.
        bonacich_coeff_var_save_eu = [bonacich_coeff_var_save_eu, nanmean(bonacich_std_eu./bonacich_mean_eu)]; % Average acrosssimulation runs s=1,...,S.
        
        %bonacich_mean_growth_eu = nanmean(bonacich_growth_eu(floor(T/2):end,:),2); % Average across simulation runs s=1,...,S.
        bonacich_std_growth_eu = std(bonacich_growth_eu(floor(T/2):end,:),0,1); % Standard deviation across time t=1,...,T.
        bonacich_std_growth_save_eu = [bonacich_std_growth_save_eu, nanmean(bonacich_std_growth_eu)]; % Average across simulation runs s=1,...,S.
 
    end
    
    %% Save statistics for different scenarios.
    n_scenario = [n_scenario n_save'];
    n_scenario_us = [n_scenario_us n_save_us'];
    n_scenario_asia = [n_scenario_asia n_save_asia'];
    n_scenario_eu = [n_scenario_eu n_save_eu'];
    bonacich_mean_scenario = [bonacich_mean_scenario bonacich_mean_save'];
    bonacich_log_sum_std_scenario = [bonacich_log_sum_std_scenario bonacich_log_sum_std_save'];
    bonacich_coeff_var_scenario = [bonacich_coeff_var_scenario bonacich_coeff_var_save'];
    bonacich_mean_scenario_us = [bonacich_mean_scenario_us bonacich_mean_save_us'];
    bonacich_log_sum_std_scenario_us = [bonacich_log_sum_std_scenario_us bonacich_log_sum_std_save_us'];
    bonacich_coeff_var_scenario_us = [bonacich_coeff_var_scenario_us bonacich_coeff_var_save_us'];
    bonacich_mean_scenario_asia = [bonacich_mean_scenario_asia bonacich_mean_save_asia'];
    bonacich_log_sum_std_scenario_asia = [bonacich_log_sum_std_scenario_asia bonacich_log_sum_std_save_asia'];
    bonacich_coeff_var_scenario_asia = [bonacich_coeff_var_scenario_asia bonacich_coeff_var_save_asia'];
    bonacich_mean_scenario_eu = [bonacich_mean_scenario_eu bonacich_mean_save_eu'];
    bonacich_log_sum_std_scenario_eu = [bonacich_log_sum_std_scenario_eu bonacich_log_sum_std_save_eu'];
    bonacich_coeff_var_scenario_eu = [bonacich_coeff_var_scenario_eu bonacich_coeff_var_save_eu'];
    out_d_std_scenario = [out_d_std_scenario out_d_std_save'];
    deg_scenario = [deg_scenario deg_save'];
    out_d_std_scenario_us = [out_d_std_scenario_us out_d_std_save_us'];
    deg_scenario_us = [deg_scenario_us deg_save_us'];
    out_d_std_scenario_asia = [out_d_std_scenario_asia out_d_std_save_asia'];
    deg_scenario_asia = [deg_scenario_asia deg_save_asia'];
    out_d_std_scenario_eu = [out_d_std_scenario_eu out_d_std_save_eu'];
    deg_scenario_eu = [deg_scenario_eu deg_save_eu'];
    if(compute_connectivity_stats)
        if(compute_weak_components) num_weak_comp_scenario = [num_weak_comp_scenario num_weak_comp_save']; end
        if(compute_weak_components) size_largest_weak_comp_scenario = [size_largest_weak_comp_scenario size_largest_weak_comp_save']; end
        if(compute_strong_components) num_strong_comp_scenario = [num_strong_comp_scenario num_strong_comp_save']; end
        if(compute_strong_components) size_largest_strong_comp_scenario = [size_largest_strong_comp_scenario size_largest_strong_comp_save']; end
        if(compute_char_path_length) char_path_length_scenario = [char_path_length_scenario char_path_length_save']; end
        if(compute_modularity) modularity_scenario = [modularity_scenario modularity_save']; end
        
        if(compute_weak_components) num_weak_comp_scenario_us = [num_weak_comp_scenario_us num_weak_comp_save_us']; end
        if(compute_weak_components) size_largest_weak_comp_scenario_us = [size_largest_weak_comp_scenario_us size_largest_weak_comp_save_us']; end
        if(compute_strong_components) num_strong_comp_scenario_us = [num_strong_comp_scenario_us num_strong_comp_save_us']; end
        if(compute_strong_components) size_largest_strong_comp_scenario_us = [size_largest_strong_comp_scenario_us size_largest_strong_comp_save_us']; end
        if(compute_char_path_length) char_path_length_scenario_us = [char_path_length_scenario_us char_path_length_save_us']; end
        if(compute_modularity) modularity_scenario_us = [modularity_scenario_us modularity_save_us']; end
        
        if(compute_weak_components) num_weak_comp_scenario_asia = [num_weak_comp_scenario_asia num_weak_comp_save_asia']; end
        if(compute_weak_components) size_largest_weak_comp_scenario_asia = [size_largest_weak_comp_scenario_asia size_largest_weak_comp_save_asia']; end
        if(compute_strong_components) num_strong_comp_scenario_asia = [num_strong_comp_scenario_asia num_strong_comp_save_asia']; end
        if(compute_strong_components) size_largest_strong_comp_scenario_asia = [size_largest_strong_comp_scenario_asia size_largest_strong_comp_save_asia']; end
        if(compute_char_path_length) char_path_length_scenario_asia = [char_path_length_scenario_asia char_path_length_save_asia']; end
        if(compute_modularity) modularity_scenario_asia = [modularity_scenario_asia modularity_save_asia']; end
        
        if(compute_weak_components) num_weak_comp_scenario_eu = [num_weak_comp_scenario_eu num_weak_comp_save_eu']; end
        if(compute_weak_components) size_largest_weak_comp_scenario_eu = [size_largest_weak_comp_scenario_eu size_largest_weak_comp_save_eu']; end
        if(compute_strong_components) num_strong_comp_scenario_eu = [num_strong_comp_scenario_eu num_strong_comp_save_eu']; end
        if(compute_strong_components) size_largest_strong_comp_scenario_eu = [size_largest_strong_comp_scenario_eu size_largest_strong_comp_save_eu']; end
        if(compute_char_path_length) char_path_length_scenario_eu = [char_path_length_scenario_eu char_path_length_save_eu']; end
        if(compute_modularity) modularity_scenario_eu = [modularity_scenario_eu modularity_save_eu']; end
    end
    
end
 
%% Plot total sales over the rewiring rate gamma.
if(multi_plot)
    subplot(3,4,1);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    %ind = find(gamma==emp_gamma);
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_mean_scenario(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ total sales [\%]');
title('All Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_average_sales_tot');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot total sales over the rewiring rate gamma for the US only.
if(multi_plot)
    subplot(3,4,2);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_mean_scenario_us(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ total sales [\%]');
title('US Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_average_sales_us');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot total sales over the rewiring rate gamma for the China only.
if(multi_plot)
    subplot(3,4,3);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_mean_scenario_asia(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ total sales [\%]');
title('Asian Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_average_sales_asia');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot total sales over the rewiring rate gamma for the China only.
if(multi_plot)
    subplot(3,4,4);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_mean_scenario_eu(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ total sales [\%]');
title('EU Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_average_sales_eu');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation in total firm sales (sum of Bonacich centralities) over the rewiring rate gamma.
if(multi_plot)
    subplot(3,4,5);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_coeff_var_scenario(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(nearest_neighbors)
        y_nn = y;
        for i=2:length(y)-1
            y_nn(i) = nanmean([y(i-1) y(i) y(i+1)]);
        end
        y = y_nn;
        hold on
    end
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ sales coeff. var. [\%]');
title('All Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_sales_tot');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation in total firm sales (sum of Bonacich centralities) over the rewiring rate gamma for the US only.
if(multi_plot)
    subplot(3,4,6);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_coeff_var_scenario_us(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(nearest_neighbors)
        ind = find(x<=gamma_plot_max);
        x = x(ind);
        y = y(ind);
        y_nn = y;
        for i=2:length(y)-1
            y_nn(i) = nanmean([y(i-1) y(i) y(i+1)]);
        end
        y = y_nn;
        hold on
    end
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ sales coeff. var. [\%]');
title('US Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_sales_us');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation in total firm sales (sum of Bonacich centralities) over the rewiring rate gamma for Asia only.
if(multi_plot)
    subplot(3,4,7);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_coeff_var_scenario_asia(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(nearest_neighbors)
        y_nn = y;
        for i=2:length(y)-1
            y_nn(i) = nanmean([y(i-1) y(i) y(i+1)]);
        end
        y = y_nn;
        hold on
    end
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ sales coeff. var. [\%]');
title('Asian Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_sales_asia');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation in total firm sales (sum of Bonacich centralities) over the rewiring rate gamma for the EU only.
if(multi_plot)
    subplot(3,4,8);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = bonacich_coeff_var_scenario_eu(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy2 = dummy1(ind);
    end
    if(~isempty(ind))
        y = 100*(dummy1-dummy2)./dummy2;
    else
        y = dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(nearest_neighbors)
        y_nn = y;
        for i=2:length(y)-1
            y_nn(i) = nanmean([y(i-1) y(i) y(i+1)]);
        end
        y = y_nn;
        hold on
    end
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ sales coeff. var. [\%]');
title('EU Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_sales_eu');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation of the out-degree over the rewiring rate gamma.
if(multi_plot)
    subplot(3,4,9);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = deg_scenario(:,c);
    dummy2 = out_d_std_scenario(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy3 = dummy1(ind);
        dummy4 = dummy2(ind);
    end
    if(~isempty(ind))
        y = 100*((dummy2./dummy1)-(dummy4./dummy3))./(dummy4./dummy3);
    else
        y = dummy2./dummy1; 
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ out-deg coeff. var. [\%]');
title('All Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_out_d_tot');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation of the out-degree over the rewiring rate gamma.
if(multi_plot)
    subplot(3,4,10);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = deg_scenario_us(:,c);
    dummy2 = out_d_std_scenario_us(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy3 = dummy1(ind);
        dummy4 = dummy2(ind);
    end
    if(~isempty(ind))
        y = 100*((dummy2./dummy1)-(dummy4./dummy3))./(dummy4./dummy3);
    else
        y = dummy2./dummy1;
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ out-deg coeff. var. [\%]');
title('US Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_out_d_us');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation of the out-degree over the rewiring rate gamma.
if(multi_plot)
    subplot(3,4,11);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = deg_scenario_asia(:,c);
    dummy2 = out_d_std_scenario_asia(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy3 = dummy1(ind);
        dummy4 = dummy2(ind);
    end
    if(~isempty(ind))
        y = 100*((dummy2./dummy1)-(dummy4./dummy3))./(dummy4./dummy3);
    else
        y = dummy2./dummy1; 
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ out-deg coeff. var. [\%]');
title('Asian Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_out_d_asia');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
%% Plot the coefficient of variation of the out-degree over the rewiring rate gamma.
if(multi_plot)
    subplot(3,4,12);
else
    figure();
end
set(0,'defaulttextinterpreter','latex')
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
for c=1:length(scenario)
    [d,ix] = min(abs(gamma-emp_gamma));
    ind = ix(1);
    dummy1 = deg_scenario_eu(:,c);
    dummy2 = out_d_std_scenario_eu(:,c);
    if(strcmp(scenario(c),'benchmark'))
        dummy3 = dummy1(ind);
        dummy4 = dummy2(ind);
    end
    if(~isempty(ind))
        y = 100*((dummy2./dummy1)-(dummy4./dummy3))./(dummy4./dummy3);
    else
        y = dummy2./dummy1; 
    end
    x = gamma;
    [x,ind] = sort(x);
    y = y(ind);
    if(strcmp(scenario(c),'benchmark'))
        plot(x,y,'ok','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--k','LineWidth', 1.5);
        hold on
    else
        plot(x,y,'or','LineWidth', 1.5)
        hold on
        y_smooth = smooth(x,y,0.95,'rloess');
        plot(x,y_smooth,'--r','LineWidth', 1.5);
        hold on
    end
end
h = vline(emp_gamma,'b','$\widehat{\gamma}$');
set(h(1),'LineWidth', 0.5)
hold on
plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
hold on
xlabel('$\gamma$');
ylabel('$\Delta$ out-deg coeff. var. [\%]');
title('EU Firms')
axis tight
xlim([gamma_plot_min gamma_plot_max]); 
hold on
set(gca,'FontSize',20);
if(~multi_plot)
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_coeff_var_out_d_eu');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
if(multi_plot)
    set(gca,'FontSize',16);
    filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_multi_plot');
    print('-depsc',strcat('Figures/',filename,'.eps'))
    print('-dpdf',strcat('Figures/',filename,'.pdf'))
end
 
if(compute_connectivity_stats)
        
    %% Plot the average out-degree over the rewiring rate gamma.
    if(multi_plot)
        figure();
        subplot(3,2,1);
    else
        figure();
    end
    set(0,'defaulttextinterpreter','latex')
    set(gca,'Layer','top')
    set(gca,'FontSize',16);
    box on
    for c=1:length(scenario)
        [d,ix] = min(abs(gamma-emp_gamma));
        ind = ix(1);
        dummy1 = deg_scenario(:,c);
        if(strcmp(scenario(c),'benchmark'))
            dummy2 = dummy1(ind);
        end
        if(~isempty(ind))
            y = 100*(dummy1-dummy2)./dummy2;
        else
            y = dummy1;
        end
        x = gamma;
        [x,ind] = sort(x);
        y = y(ind);
        if(strcmp(scenario(c),'benchmark'))
            plot(x,y,'ok','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--k','LineWidth', 1.5);
            hold on
        else
            plot(x,y,'or','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--r','LineWidth', 1.5);
            hold on
        end
    end
    h = vline(emp_gamma,'b','$\widehat{\gamma}$');
    set(h(1),'LineWidth', 0.5)
    hold on
    plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
    hold on
    xlabel('$\gamma$');
    ylabel('$\Delta$ av. deg. [\%]');
    title('All Firms')
    axis tight
    xlim([gamma_plot_min gamma_plot_max]); 
    hold on
    set(gca,'FontSize',16);
    if(~multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_av_out_d');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,2);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_weak_comp_scenario(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. weak comp. [\%]');
        %ylabel('num. weak comp.');
        title('All Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_weak_comp');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,3);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_weak_comp_scenario(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest weak comp. [\%]');
        title('All Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_weak_comp');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,4);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_strong_comp_scenario(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. strong comp. [\%]');
        title('All Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_strong_comp');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,5);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_strong_comp_scenario(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest strong comp. [\%]');
        title('All Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_strong_comp');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_char_path_length) 
        if(multi_plot)
            subplot(3,2,6);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = char_path_length_scenario(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ char. path length [\%]');
        title('All Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_char_path_length');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_modularity) % modularity_scenario
        if(multi_plot)
            subplot(4,2,7);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = modularity_scenario(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ modularity [\%]');
        %ylabel('modularity');
        title('All Firms')
axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_modularity');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_multi_plot_connectivity_stats');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
end
 
if(compute_connectivity_stats)
        
    %% Plot the average out-degree over the rewiring rate gamma.
    if(multi_plot)
        figure();
        subplot(3,2,1);
    else
        figure();
    end
    set(0,'defaulttextinterpreter','latex')
    set(gca,'Layer','top')
    set(gca,'FontSize',16);
    box on
    for c=1:length(scenario)
        [d,ix] = min(abs(gamma-emp_gamma));
        ind = ix(1);
        dummy1 = deg_scenario_us(:,c);
        if(strcmp(scenario(c),'benchmark'))
            dummy2 = dummy1(ind);
        end
        if(~isempty(ind))
            y = 100*(dummy1-dummy2)./dummy2;
        else
            y = dummy1;
        end
        x = gamma;
        [x,ind] = sort(x);
        y = y(ind);
        if(strcmp(scenario(c),'benchmark'))
            plot(x,y,'ok','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--k','LineWidth', 1.5);
            hold on
        else
            plot(x,y,'or','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--r','LineWidth', 1.5);
            hold on
        end
    end
    h = vline(emp_gamma,'b','$\widehat{\gamma}$');
    set(h(1),'LineWidth', 0.5)
    hold on
    plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
    hold on
    xlabel('$\gamma$');
    ylabel('$\Delta$ av. deg. [\%]');
    title('US Firms')
    axis tight
    xlim([gamma_plot_min gamma_plot_max]); 
    hold on
    set(gca,'FontSize',16);
    if(~multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_av_out_d_us');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,2);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_weak_comp_scenario_us(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. weak comp. [\%]');
        title('US Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_weak_comp_us');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,3);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_weak_comp_scenario_us(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest weak comp. [\%]');
        title('US Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_weak_comp_us');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,4);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_strong_comp_scenario_us(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. strong comp. [\%]');
        title('US Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_strong_comp_us');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,5);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_strong_comp_scenario_us(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest strong comp. [\%]');
        title('US Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_strong_comp_us');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_char_path_length) 
        if(multi_plot)
            subplot(3,2,6);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = char_path_length_scenario_us(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ char. path length [\%]');
        title('US Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_char_path_length_us');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_modularity) % modularity_scenario_us
        if(multi_plot)
            subplot(4,2,7);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            %ind = find(gamma==emp_gamma);
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = modularity_scenario_us(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ modularity [\%]');
        %ylabel('modularity');
        title('US Firms')
axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_modularity_us');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_multi_plot_connectivity_stats_us');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
end
 
if(compute_connectivity_stats)
        
    %% Plot the average out-degree over the rewiring rate gamma.
    if(multi_plot)
        figure();
        subplot(3,2,1);
    else
        figure();
    end
    set(0,'defaulttextinterpreter','latex')
    set(gca,'Layer','top')
    set(gca,'FontSize',16);
    box on
    for c=1:length(scenario)
        [d,ix] = min(abs(gamma-emp_gamma));
        ind = ix(1);
        dummy1 = deg_scenario_asia(:,c);
        if(strcmp(scenario(c),'benchmark'))
            dummy2 = dummy1(ind);
        end
        if(~isempty(ind))
            y = 100*(dummy1-dummy2)./dummy2;
        else
            y = dummy1;
        end
        x = gamma;
        [x,ind] = sort(x);
        y = y(ind);
        if(strcmp(scenario(c),'benchmark'))
            plot(x,y,'ok','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--k','LineWidth', 1.5);
            hold on
        else
            plot(x,y,'or','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--r','LineWidth', 1.5);
            hold on
        end
    end
    h = vline(emp_gamma,'b','$\widehat{\gamma}$');
    set(h(1),'LineWidth', 0.5)
    hold on
    plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
    hold on
    xlabel('$\gamma$');
    ylabel('$\Delta$ av. deg. [\%]');
    %ylabel('av. deg.');
    title('Asian Firms')
    axis tight
    xlim([gamma_plot_min gamma_plot_max]); 
    hold on
    set(gca,'FontSize',16);
    if(~multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_av_out_d_asia');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,2);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_weak_comp_scenario_asia(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. weak comp. [\%]');
        title('Asian Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_weak_comp_asia');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_weak_components)
        if(multi_plot)
            subplot(3,2,3);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_weak_comp_scenario_asia(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest weak comp. [\%]');
        title('Asian Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_weak_comp_asia');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,4);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_strong_comp_scenario_asia(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(nearest_neighbors)
                y_nn = y;
                for i=2:length(y)-1
                    y_nn(i) = nanmean([y(i-1) y(i) y(i+1)]);
                end
                y = y_nn;
                hold on
            end
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. strong comp. [\%]');
        title('Asian Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_strong_comp_asia');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,5);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_strong_comp_scenario_asia(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest strong comp. [\%]');
        title('Asian Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_strong_comp_asia');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_char_path_length) 
        if(multi_plot)
            subplot(3,2,6);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = char_path_length_scenario_asia(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ char. path length [\%]');
        title('Asian Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_char_path_length_asia');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_modularity) % modularity_scenario_asia
        if(multi_plot)
            subplot(4,2,7);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            %ind = find(gamma==emp_gamma);
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = modularity_scenario_asia(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ modularity [\%]');
        %ylabel('modularity');
        title('Asian Firms')
axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_modularity_asia');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_multi_plot_connectivity_stats_asia');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
end
 
if(compute_connectivity_stats)
        
    %% Plot the average out-degree over the rewiring rate gamma.
    if(multi_plot)
        figure();
        subplot(3,2,1);
    else
        figure();
    end
    set(0,'defaulttextinterpreter','latex')
    set(gca,'Layer','top')
    set(gca,'FontSize',16);
    box on
    for c=1:length(scenario)
        [d,ix] = min(abs(gamma-emp_gamma));
        ind = ix(1);
        dummy1 = deg_scenario_eu(:,c);
        if(strcmp(scenario(c),'benchmark'))
            dummy2 = dummy1(ind);
        end
        if(~isempty(ind))
            y = 100*(dummy1-dummy2)./dummy2;
        else
            y = dummy1;
        end
        x = gamma;
        [x,ind] = sort(x);
        y = y(ind);
        if(strcmp(scenario(c),'benchmark'))
            plot(x,y,'ok','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--k','LineWidth', 1.5);
            hold on
        else
            plot(x,y,'or','LineWidth', 1.5)
            hold on
            y_smooth = smooth(x,y,0.95,'rloess');
            plot(x,y_smooth,'--r','LineWidth', 1.5);
            hold on
        end
    end
    h = vline(emp_gamma,'b','$\widehat{\gamma}$');
    set(h(1),'LineWidth', 0.5)
    hold on
    plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
    hold on
    xlabel('$\gamma$');
    ylabel('$\Delta$ av. deg. [\%]');
    title('EU Firms')
    axis tight
    xlim([gamma_plot_min gamma_plot_max]); 
    hold on
    set(gca,'FontSize',16);
    if(~multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_av_out_d_eu');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,2);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_weak_comp_scenario_eu(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. weak comp. [\%]');
        title('EU Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_weak_comp_eu');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_weak_components) 
        if(multi_plot)
            subplot(3,2,3);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_weak_comp_scenario_eu(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest weak comp. [\%]');
        title('EU Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_weak_comp_eu');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,4);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = num_strong_comp_scenario_eu(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ num. strong comp. [\%]');
        title('EU Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_num_strong_comp_eu');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_strong_components) 
        if(multi_plot)
            subplot(3,2,5);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = size_largest_strong_comp_scenario_eu(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ size largest strong comp. [\%]');
        title('EU Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_size_largest_strong_comp_eu');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_char_path_length) 
        if(multi_plot)
            subplot(3,2,6);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = char_path_length_scenario_eu(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ char. path length [\%]');
        title('EU Firms')
        axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_char_path_length_eu');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(compute_modularity) % modularity_scenario_eu
        if(multi_plot)
            subplot(4,2,7);
        else
            figure();
        end
        set(0,'defaulttextinterpreter','latex')
        set(gca,'Layer','top')
        set(gca,'FontSize',16);
        box on
        for c=1:length(scenario)
            %ind = find(gamma==emp_gamma);
            [d,ix] = min(abs(gamma-emp_gamma));
            ind = ix(1);
            dummy1 = modularity_scenario_eu(:,c);
            if(strcmp(scenario(c),'benchmark'))
                dummy2 = dummy1(ind);
            end
            if(~isempty(ind))
                y = 100*(dummy1-dummy2)./dummy2;
            else
                y = dummy1;
            end
            x = gamma;
            [x,ind] = sort(x);
            y = y(ind);
            if(strcmp(scenario(c),'benchmark'))
                plot(x,y,'ok','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--k','LineWidth', 1.5);
                hold on
            else
                plot(x,y,'or','LineWidth', 1.5)
                hold on
                y_smooth = smooth(x,y,0.95,'rloess');
                plot(x,y_smooth,'--r','LineWidth', 1.5);
                hold on
            end
        end
        h = vline(emp_gamma,'b','$\widehat{\gamma}$');
        set(h(1),'LineWidth', 0.5)
        hold on
        plot(x,zeros(length(x),1),'-b','LineWidth', 0.5)
        hold on
        xlabel('$\gamma$');
        ylabel('$\Delta$ modularity [\%]');
        %ylabel('modularity');
        title('EU Firms')
axis tight
        xlim([gamma_plot_min gamma_plot_max]); 
        hold on
        set(gca,'FontSize',16);
        if(~multi_plot)
            filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_modularity_eu');
            print('-depsc',strcat('Figures/',filename,'.eps'))
            print('-dpdf',strcat('Figures/',filename,'.pdf'))
        end
    end
    if(multi_plot)
        set(gca,'FontSize',16);
        filename = strcat('Prodnet','_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(mean(delta)),'_zeta_',num2str(mean(zeta)),'_min_gamma_',num2str(min(gamma)),'_max_gamma_',num2str(max(gamma)),'_rho_',num2str(mean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_autarky_gamma_multi_plot_connectivity_stats_eu');
        print('-depsc',strcat('Figures/',filename,'.eps'))
        print('-dpdf',strcat('Figures/',filename,'.pdf'))
    end
end
 
t = toc;
disp(['Elapsed time is ' datestr(datenum(0,0,0,0,0,t),'HH:MM:SS')])
